
function sandbox = fsb_calculate_red_map(sandbox,idat,statsel,slice_n,hrf_pred,trial)

% FSB : Calculate statistical map for displayed slices only
%
% EXAMPLE:
% sandbox = fsb_calculate_red_map(sandbox,idat,'corrc'l,10,10,10,1)
%
% INPUT:
% sandbox:      sandbox experiment struct
% idat:             4-D image data
% statsel:         Statistical model
% slice_n : Selected slices
% hrf_pred:     Selected hemodynamic predictor
% trial: vector with trial inclusion information (logical,can be omitted)
%
% OUTPUT:
% sandbox: annotated sandbox struct with activation maps for the slices selected
%
% CALLED BY:
% FSB.m
%
% NOTES:
% Calculates a statistical map for the displayed sections/slices only
% Used for fast map display in fMRI Sandbox
%
% Copyright 2010 MPI for Biological Cybernetics
% Author: Steffen Stoewer
% License:GNU GPL, no express or implied warranties
%
% $ Revision 1.0
%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

[a b c d] = size(idat);
sandbox.hemodynamics = single(sandbox.hemodynamics);
spare_hemodynamics = sandbox.hemodynamics;

idat = single(idat);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Check if trial information provided, and correct hemodynamics information
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

if nargin == 6
    sandbox.hemodynamics = sandbox.hemodynamics(trial,:);
end

switch statsel

    case 'glm';
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        % Do calculation for General Linear Model
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        h = waitbar(0,'Calculating reduced GLM map...');
        for x =1:a;
            waitbar((x/3)/a)
            for y=1:b;
                if sandbox.stats.glm_prob(x,y,slice_n(3),hrf_pred)==0;
                    [besti2,dev,statsi] = glmfit(sandbox.hemodynamics,squeeze(idat(x,y,slice_n(3),:)));
                    sandbox.stats.glm_prob(x,y,slice_n(3),:) = statsi.p(2:end);
                end
            end

            for z = 1:c
                if sandbox.stats.glm_prob(x,slice_n(2),z,hrf_pred)==0;
                    [besti2,dev,statsi] = glmfit(sandbox.hemodynamics,squeeze(idat(x,slice_n(2),z,:)));
                    sandbox.stats.glm_prob(x,slice_n(2),z,:) = statsi.p(2:end);
                end
            end
        end

        for y=1:b;
            waitbar(((y/3)+(2*b)/3)/b)
            for z=1:c;
                if sandbox.stats.glm_prob(slice_n(1),y,z,hrf_pred)==0;
                    [besti2,dev,statsi] = glmfit(sandbox.hemodynamics,squeeze(idat(slice_n(1),y,z,:)));
                    sandbox.stats.glm_prob(slice_n(1),y,z,:) = statsi.p(2:end);
                end
            end
        end
        close(h);

    case 'corrc';
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        % Do calculation for Correlation Coefficients
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

        try % This is the faster codepath, but does not show a waitbar

            sandbox.hrf_pred = hrf_pred;
            sandbox.stats.r_cor = fsb_calculate_corr(idat,sandbox,slice_n);

        catch
            h = waitbar(0,'Calculating reduced correlation map...');
            for x =1:a;
                waitbar((x/3)/a)
                for y=1:b;
                    if sandbox.stats.r_cor(x,y,slice_n(3),hrf_pred)==0;
                        R= corrcoef(squeeze(idat(x,y,slice_n(3),:)),sum(sandbox.hemodynamics(:,hrf_pred),2));
                        sandbox.stats.r_cor(x,y,slice_n(3),hrf_pred) = R(2,1);
                    end
                end

                for z = 1:c
                    if sandbox.stats.r_cor(x,slice_n(2),z,hrf_pred)==0;
                        R = corrcoef(squeeze(idat(x,slice_n(2),z,:)),sum(sandbox.hemodynamics(:,hrf_pred),2));
                        sandbox.stats.r_cor(x,slice_n(2),z,hrf_pred) = R(2,1);
                    end
                end
            end


            for y=1:b;
                waitbar(((y/3)+(2*b)/3)/b)
                for z=1:c;
                    if sandbox.stats.r_cor(slice_n(1),y,z,hrf_pred)==0;
                        R= corrcoef(squeeze(idat(slice_n(1),y,z,:)),sum(sandbox.hemodynamics(:,hrf_pred),2));
                        sandbox.stats.r_cor(slice_n(1),y,z,hrf_pred) = R(2,1);
                    end
                end
            end
            close(h)
        end


    case 'robreg';
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        % Do calculation for robust regression
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        h = waitbar(0,'Calculating reduced robust regression map...');
        for x =1:a;
            waitbar((x/3)/a)
            for y=1:b;
                if sandbox.stats.rob_prob(x,y,slice_n(3),hrf_pred)==0;
                    [robr,statsi] = robustfit(sandbox.hemodynamics,squeeze(idat(x,y,slice_n(3),:)));
                    sandbox.stats.rob_est(x,y,slice_n(3),:)=robr(2:end);
                    sandbox.stats.rob_prob(x,y,slice_n(3),:) = statsi.p(2:end);
                end
            end

            for z = 1:c
                if sandbox.stats.rob_prob(x,slice_n(2),z,hrf_pred)==0;
                    [robr,statsi] = robustfit(sandbox.hemodynamics,squeeze(idat(x,slice_n(2),z,:)));
                    sandbox.stats.rob_est(x,slice_n(2),z,:)=robr(2:end);
                    sandbox.stats.rob_prob(x,slice_n(2),z,:) = statsi.p(2:end);
                end
            end
        end

        for y=1:b;
            waitbar(((y/3)+(2*b)/3)/b)
            for z=1:c;
                if sandbox.stats.rob_prob(slice_n(1),y,z,hrf_pred)==0;
                    [robr,statsi] = robustfit(sandbox.hemodynamics,squeeze(idat(slice_n(1),y,z,:)));
                    sandbox.stats.rob_est(slice_n(1),y,z,:)=robr(2:end);
                    sandbox.stats.rob_prob(slice_n(1),y,z,:) = statsi.p(2:end);
                end
            end
        end
        close(h);
end

sandbox.hemodynamics = spare_hemodynamics;

end
